https://rspatial.org/raster/rs/2-exploration.html
unzip('data/rsdata.zip', exdir='data')
## Warning in unzip("data/rsdata.zip", exdir = "data"): error 1 in extracting from
## zip file
“Create RasterLayer objects for single Landsat layers (bands)”
library(raster)
## Loading required package: sp
# Blue
b2 <- raster('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- raster('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- raster('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- raster('data/rs/LC08_044034_20170614_B5.tif')
see the spatial resolution, extent, number of layers, coordinate reference system, etc
b2
## class : RasterLayer
## dimensions : 1245, 1497, 1863765 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : LC08_044034_20170614_B2.tif
## names : LC08_044034_20170614_B2
## values : 0.0748399, 0.7177562 (min, max)
b3
## class : RasterLayer
## dimensions : 1245, 1497, 1863765 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : LC08_044034_20170614_B3.tif
## names : LC08_044034_20170614_B3
## values : 0.04259216, 0.6924697 (min, max)
b4
## class : RasterLayer
## dimensions : 1245, 1497, 1863765 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : LC08_044034_20170614_B4.tif
## names : LC08_044034_20170614_B4
## values : 0.02084067, 0.7861769 (min, max)
b5
## class : RasterLayer
## dimensions : 1245, 1497, 1863765 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : LC08_044034_20170614_B5.tif
## names : LC08_044034_20170614_B5
## values : 0.0008457669, 1.012432 (min, max)
# coordinate reference system (CRS)
crs(b2)
## CRS arguments:
## +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b3)
## CRS arguments:
## +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b4)
## CRS arguments:
## +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b5)
## CRS arguments:
## +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
# Number of cells, rows, columns
ncell(b2)
## [1] 1863765
dim(b2)
## [1] 1245 1497 1
# spatial resolution
res(b2)
## [1] 30 30
# Number of bands
nlayers(b2)
## [1] 1
# Do the bands have the same extent, number of rows and columns, projection, resolution, and origin
compareRaster(b2,b3)
## [1] TRUE
s <- stack(b5, b4, b3)
s
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B5, LC08_044034_20170614_B4, LC08_044034_20170614_B3
## min values : 0.0008457669, 0.0208406653, 0.0425921641
## max values : 1.0124315, 0.7861769, 0.6924697
plot(s)
plot(b5)
# create list of raster objects
filenames <- paste0("data/rs/LC08_044034_20170614_B", 1:11, ".tif")
filenames
## [1] "data/rs/LC08_044034_20170614_B1.tif"
## [2] "data/rs/LC08_044034_20170614_B2.tif"
## [3] "data/rs/LC08_044034_20170614_B3.tif"
## [4] "data/rs/LC08_044034_20170614_B4.tif"
## [5] "data/rs/LC08_044034_20170614_B5.tif"
## [6] "data/rs/LC08_044034_20170614_B6.tif"
## [7] "data/rs/LC08_044034_20170614_B7.tif"
## [8] "data/rs/LC08_044034_20170614_B8.tif"
## [9] "data/rs/LC08_044034_20170614_B9.tif"
## [10] "data/rs/LC08_044034_20170614_B10.tif"
## [11] "data/rs/LC08_044034_20170614_B11.tif"
“The layers represent reflection intensity in the following wavelengths: Ultra Blue, Blue, Green, Red, Near Infrared (NIR), Shortwave Infrared (SWIR) 1, Shortwave Infrared (SWIR) 2, Panchromatic, Cirrus, Thermal Infrared (TIRS) 1, Thermal Infrared (TIRS) 2.
We won’t use the last four layers and you will see how to remove those in following sections.”
landsat <- stack(filenames)
landsat
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 11 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11
## min values : 9.641791e-02, 7.483990e-02, 4.259216e-02, 2.084067e-02, 8.457669e-04, -7.872183e-03, -5.052945e-03, 3.931751e-02, -4.337332e-04, 2.897978e+02, 2.885000e+02
## max values : 0.73462820, 0.71775615, 0.69246972, 0.78617686, 1.01243150, 1.04320455, 1.11793602, 0.82673049, 0.03547901, 322.43139648, 317.99530029
plot individual layers in a 2x2 grid
legend values between 0-1
“different surface features reflect the incident solar radiation differently.
Each layer represent how much incident solar radiation is reflected for a particular wavelength range.
For example, vegetation reflects more energy in NIR than other wavelengths and thus appears brighter.
In contrast, water absorbs most of the energy in the NIR wavelength and it appears dark.”
par(mfrow = c(2, 2))
plot(b2, main = "Blue", col = gray(0:100 / 100))
plot(b3, main = "Green", col = gray(0:100 / 100))
plot(b4, main = "Red", col = gray(0:100 / 100))
plot(b5, main = "NIR", col = gray(0:100 / 100))
“(vegetation in green, water blue etc), we need bands in the red, green and blue regions.
For this Landsat image, band 4 (red), 3 (green), and 2 (blue) can be used.”
par(oma = c(2, 0, 1, 4)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.8, 3.1, .55, 0.35))
landsatRGB <- stack(b4, b3, b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
## False colour image
par(oma = c(2, 0, 1, 4)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.8, 3.1, .55, 0.35))
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
par(mfrow = c(1, 2))
par(oma = c(0, 0.5, 0, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "s") # "s" square plotting region vs m - maximal
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
help(plotRGB)
# select first 3 bands
landsat
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 11 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11
## min values : 9.641791e-02, 7.483990e-02, 4.259216e-02, 2.084067e-02, 8.457669e-04, -7.872183e-03, -5.052945e-03, 3.931751e-02, -4.337332e-04, 2.897978e+02, 2.885000e+02
## max values : 0.73462820, 0.71775615, 0.69246972, 0.78617686, 1.01243150, 1.04320455, 1.11793602, 0.82673049, 0.03547901, 322.43139648, 317.99530029
landsatsub1 <- subset(landsat, 4:2)
landsatsub1
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2
## min values : 0.02084067, 0.04259216, 0.07483990
## max values : 0.7861769, 0.6924697, 0.7177562
# same as
landsatsub2 <- landsat[[4:2]]
landsatsub2
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2
## min values : 0.02084067, 0.04259216, 0.07483990
## max values : 0.7861769, 0.6924697, 0.7177562
# Number of bands in the original and new data
nlayers(landsat)
## [1] 11
nlayers(landsatsub1)
## [1] 3
nlayers(landsatsub2)
## [1] 3
landsat <- subset(landsat, 1:7)
landsat
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7
## min values : 0.0964179114, 0.0748399049, 0.0425921641, 0.0208406653, 0.0008457669, -0.0078721829, -0.0050529451
## max values : 0.7346282, 0.7177562, 0.6924697, 0.7861769, 1.0124315, 1.0432045, 1.1179360
names(landsat)
## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
names(landsat)
## [1] "ultra.blue" "blue" "green" "red" "NIR"
## [6] "SWIR1" "SWIR2"
“used to limit analysis to a geographic subset of the image”
extent(landsat)
## class : Extent
## xmin : 594090
## xmax : 639000
## ymin : 4190190
## ymax : 4227540
e <- extent(624387, 635752, 4200047, 4210939)
e
## class : Extent
## xmin : 624387
## xmax : 635752
## ymin : 4200047
## ymax : 4210939
# crop landsat by the extent
landsatcrop <- crop(landsat, e)
landsatcrop
## class : RasterBrick
## dimensions : 363, 379, 137577, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 624390, 635760, 4200060, 4210950 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : memory
## names : ultra.blue, blue, green, red, NIR, SWIR1, SWIR2
## min values : 0.101796150, 0.075989284, 0.044045158, 0.030556191, 0.007243267, 0.001865030, 0.003144530
## max values : 0.4548080, 0.4746728, 0.5034941, 0.5592065, 0.7013395, 0.9038041, 0.9750224
par(mfrow = c(1, 2))
landsatsub3 <- subset(landsat, 4:2)
landsatsub4 <- subset(landsat, 5:3)
plotRGB(landsatsub3, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
plotRGB(landsatsub4, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
par(ps = 9) # point text size
par(oma = c(0, 1.5, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(1.2, 4.2, .55, 0.35))
par(mar = c(5, 5, 4, 2))
# ?par # help for par() function
landsatsub5 <- subset(landsatcrop, 4:2)
plotRGB(landsatsub5, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nTrue Colour Composite")
par(ps = 9) # point text size
par(oma = c(0, 1.5, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(1.2, 4.2, .55, 0.35))
par(mar = c(5, 5, 4, 2))
# ?par # help for par() function
# title(sub = "Comparing True and False Colour Landsat Composites") # only works for area of one graph
# xlab = "UTM E/W", ylab = "UTM N/S", takes too much space
landsatsub6 <- subset(landsatcrop, 5:3)
plotRGB(landsatsub6, axes = TRUE, stretch = "lin", main = "Landsat\nFalse Colour Composite")
par(mfrow = c(1, 2))
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
# ?par # help for par() function
landsatsub5 <- subset(landsatcrop, 4:2)
plotRGB(landsatsub5, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nTrue Colour Composite")
# title(sub = "Comparing True and False Colour Landsat Composites") # only works for area of one graph
landsatsub6 <- subset(landsatcrop, 5:3)
plotRGB(landsatsub6, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nFalse Colour Composite")
jpeg(file="landsatRGB-TCC.jpeg")
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
dev.off()
## quartz_off_screen
## 2
jpeg(file="landsatRGB-FCC.jpeg")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
dev.off()
## quartz_off_screen
## 2
tiff(file="landsatRGB-TCC.tif")
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
dev.off()
## quartz_off_screen
## 2
tiff(file="landsatRGB-FCC.tif")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
dev.off()
## quartz_off_screen
## 2
rf <- writeRaster(landsatcrop, filename = "cropped-landsat.tif", overwrite=TRUE)
rf <- writeRaster(landsatsub5, filename = "cropped-landsat5.tif", overwrite=TRUE)
rf <- writeRaster(s, filename = "cropped-landsat-TCC.tif", overwrite=TRUE)
# S4 method for RasterStackBrick,character
rf <- writeRaster(landsatFCC, filename=file.path("landsatFCC.tif"), format="GTiff", overwrite=TRUE)
rf <- writeRaster(landsatRGB, filename=file.path("landsatTCC.tif"), format="GTiff", overwrite=TRUE)
class(landsatcrop)
## [1] "RasterBrick"
## attr(,"package")
## [1] "raster"
class(landsat)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"
There is a high correlation (99%) between the ultra-blue and blue wavelengths.
so not much info will be lost if sing only one of these bands
pairs(landsatcrop[[1:2]], main = "Correlation between Ultra-blue versus Blue")
A triangle distribution. Vegetation reflects higher in NIR than red wavelengths
“and creates the upper corner close to NIR (y) axis.
Water absorbs energy from all the bands and occupies the location close to origin.
The furthest corner is created due to highly reflecting surface features like bright soil or concrete.”
pairs(landsatcrop[[4:5]], main = "Red vs NIR")
# load the polygons with land use land cover information
samp <- readRDS('data/rs/samples.rds')
# generate 300 point samples from the polygons
ptsamp <- spsample(samp, 300, type='regular')
## Warning in proj4string(obj): CRS object has comment, which is lost in output
# add the land cover class to the points
ptsamp$class <- over(ptsamp, samp)$class
# extract values with points
df <- extract(landsat, ptsamp)
# To see some of the reflectance values
head(df)
## ultra.blue blue green red NIR SWIR1 SWIR2
## [1,] 0.1294681 0.1235693 0.12922950 0.17581198 0.3241255 0.3347085 0.1969346
## [2,] 0.1422847 0.1373185 0.14692563 0.18604797 0.2774130 0.3170558 0.2110958
## [3,] 0.1327861 0.1148514 0.09830463 0.09570226 0.1651856 0.2010984 0.1640145
## [4,] 0.1344559 0.1170851 0.10006124 0.09889016 0.1687422 0.2109006 0.1755517
## [5,] 0.1311379 0.1291861 0.14445338 0.19541651 0.3324314 0.3122631 0.1759638
## [6,] 0.1384245 0.1338053 0.14467025 0.19270571 0.3112655 0.3450096 0.2207463
ms <- aggregate(df, list(ptsamp$class), mean)
# instead of the first column, we use row names
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms
## ultra.blue blue green red NIR SWIR1
## built 0.1749001 0.16639696 0.16668714 0.17878095 0.22796993 0.23626448
## cropland 0.1118587 0.08961921 0.08389110 0.05294961 0.46342623 0.15512600
## fallow 0.1321714 0.11686468 0.10508105 0.11573817 0.18102815 0.23419503
## open 0.1385035 0.13693335 0.15151051 0.20491259 0.34288999 0.35559024
## water 0.1337126 0.11681711 0.09981322 0.07909337 0.04886901 0.03346269
## SWIR2
## built 0.19687054
## cropland 0.06828988
## fallow 0.19450095
## open 0.21181017
## water 0.02714454
"The spectral profile plots reflectance of different earth’s surface features
‘Water’ has low reflection in all wavelengths
‘built’, ‘fallow’ and ‘open’ have have high reflectance in longer wavelengths
# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")
# add the different classes
for (i in 1:nrow(ms)){
lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Title
title(main="Spectral Profile from Landsat", font.main = 2)
# Legend
legend("topleft", rownames(ms),
cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")
rasterperformed on each pixel / grid cell
raslist <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <-stack(raslist)
landsatRGB <- landsat[[c(4, 3, 2)]]
landsatFCC <- landsat[[c(5, 4, 3)]]
vi <- function(img, k, i) {
bk <- img[[k]]
bi <- img[[i]]
vi <- (bk - bi) / (bk + bi)
return(vi)
}
to view greenness
Landsat bands NIR = 5, red = 4 #### NDVI one way
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Landsat NDVI", sub = "Level of greenness")
vi2 <- function(x, y){
(x - y) / (x + y)
}
ndvi2 <- overlay(landsat[[5]], landsat[[4]], fun = vi2)
plot(ndvi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Landsat NDVI", sub = "Level of greenness")
see spectral profile for bands with max / min reflectance
vi <- function(img, k, i) {
bk <- img[[k]]
bi <- img[[i]]
vi <- (bk - bi) / (bk + bi)
return(vi)
}
Landsat bands: green = 3 and NIR = 5
ndwi <- vi(landsat, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Body Index: Landsat Bands green (3) and NIR (5)", sub = "Monitor water in water bodies")
Landsat bands: NIR = 5 and SWIR = 7
(Note: very similar to Built Environment Index 2 Bands: NIR = 5 and green = 3)
ndwi2 <- vi(landsat, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (7)", sub = "Monitor change in leaf water content")
Landsat bands: NIR = 5 and SWIR = 6
this is not the band to identify water in leaves
ndwi2 <- vi(landsat, 5, 6)
plot(ndwi2, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (6)", sub = "Monitor change in leaf water content")
see spectral profile for bands with max / min reflectance
e.g. Identify built areas NDBI (Normalized Difference Built-up Index)
Landsat bands: 4 = red, 2 = green, 7 = SWIR
https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
bk <- img[[k]]
bi <- img[[i]]
bp <- img[[p]]
bi <- (bk - bi) * bp / (bk + bi) * bp
return(bi)
}
Landsat Bands: 4 = red, 2 = blue, 7 = SWIR
https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm
ndbi <- bi(landsat, 4, 2, 7)
plot(ndbi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands red (4), blue (2), and SWIR (7)", sub = "Monitor change in built environment (4-2)*7/(4+2)*7")
Landsat bands: 6 = SWIR, 5 = NIR
https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i) {
bk <- img[[k]]
bi <- img[[i]]
bi <- (bk - bi) / (bk + bi)
return(bi)
}
ndbi <- bi(landsat, 6, 5)
plot(ndbi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands SWIR (6) and NIR (5)", sub = "Monitor change in built environment (6-5)/(6+5)")
Landsat bands: 4 = red, 5 = NIR
** Don’t Use (Note: too similar to Built Environment Index 2 Bands: NIR = 5 and green = 3) https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm
bei <- vi(landsat, 4, 5)
plot(bei, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands NIR (5) and green (3)", sub = "Monitor change in built environment")
par(mfrow = c(2, 2)) # rows, columns
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.35, 0.35, 0.35, 0.35))
# ?par # help for par() function
###############
# NDVI plot
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat NDVI Greeness Index")
###############
# NDWI plot leaf water
ndwi2 <- vi(landsat, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (7)")
###############
# NDWI plot water bodies
ndwi <- vi(landsat, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), main = "Water Body Index: Landsat Bands green (3) and NIR (5)")
###############
# Built environment index
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
bk <- img[[k]]
bi <- img[[i]]
bp <- img[[p]]
bi <- (bk - bi) * bp / (bk + bi) * bp
return(bi)
}
bsi <- bi(landsat, 4, 2, 7)
plot(bsi, col = rev(terrain.colors(10)), main = "Built Index: Landsat Bands Red, Green and SWIR [(4-2)*7/(4+2)*7]")
par(mfrow = c(2, 2)) # rows, columns
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.35, 0.35, 0.35, 0.35))
# ?par # help for par() function
###############
# NDVI plot
ndvi <- vi(landsatcrop, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat NDVI Greeness Index")
###############
# NDWI plot leaf water
ndwi2 <- vi(landsatcrop, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (7)")
###############
# NDWI plot water bodies
ndwi <- vi(landsatcrop, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), main = "Water Body Index: Landsat Bands green (3) and NIR (5)")
###############
# Built environment index
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
bk <- img[[k]]
bi <- img[[i]]
bp <- img[[p]]
bi <- (bk - bi) * bp / (bk + bi) * bp
return(bi)
}
bsi <- bi(landsatcrop, 4, 2, 7)
plot(bsi, col = rev(terrain.colors(10)), main = "Built Index: Landsat Bands Red, Green and SWIR [(4-2)*7/(4+2)*7]")
###############
# Plot Spectra
# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")
# add the different classes
for (i in 1:nrow(ms)){
lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Title
title(main="Spectral Profile from Landsat", font.main = 2)
# Legend
legend("topleft", rownames(ms),
cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")